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ABSTRACT 

We use a probabilistic method to compute the propagation of an ionization front 
corresponding to the re-ionization of the intergalactic medium in a ACDM cosmology, 
including both hydrogen and helium. The effects of radiative transfer substantially 
boost the temperature of the ionized gas over the case of uniform re-ionization. The 
resulting temperature-density relation of the ionized gas is both non-monotonic and 
multiple- valued, reflecting the non-local character of radiative transfer and suggesting 
that a single polytropic relation between local gas density and temperatue is a poor 
description of the thermodynamic state of baryons in the post-reionization universe. 
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1 INTRODUCTION 

Hydrodynamical simulations of structure formation in the 
universe have led to fundamental insights into the struc- 
ture and evolution of the Intergalactic Medium (IGM) (Cen 
et al. 1994; Zhang, Anninos & Norman 1995; Hernquist 
et al. 1996; Zhang et al. 1997; Bond & Wadsley 1997; 
Theuns, Leonard & Efstathiou 1998). Comparisons with the 
Lya forest, as measured in high redshift Quasi-Stellar Ob- 
ject (QSO) spectra, show that the simulations broadly repro- 
duce the statistical properties of the IGM. Precision compar- 
isons with the highest spectral resolution measurements re- 
cover the cumulative flux distributions and H i column den- 
sity distributions to an accuracy of a few percent (Meiksin, 
Bryan & Machacek 2001). More problematic are the pre- 
dicted widths of the absorption features, which appear 
to require additional sources of broadening, perhaps late 
He II re-ionization, to match the measured widths (Theuns 
et al. 1999; Bryan & Machacek 2000; Meiksin et al. 2001). 

Many of the properties of the IGM may be attributed to 
the gravitational instability of the dark matter alone, as the 
baryon density fluctuations closely follow those of the dark 
matter (Zhang et al. 1998). This has led to the development 
of pseudo-hydrodynamic simulations based on pure gravity 
in which the baryon fluctuations are assumed to exactly fol- 
low the dark matter and the temperature to be derived from 
an effective equation of state (Petitjean, Miicket & Kates 
1995; Croft et al. 1998; Gnedin & Hui 1998; Meiksin & 
White 2001). In recent years, such simulations have been in- 
creasingly relied on for predicting the flux power spectrum 
of the Lya forest (Meiksin & White 2001; Zaldarriaga, Hui 
& Tegmark 2001; Croft et al. 2002; Meiksin & White 2003b; 
Seljak, McDonald & Makarov 2003). 



The simulations have generally been done assum- 
ing sudden early homogeneous H i and He ii re-ionization. 
While tentative steps have been taken to introduce ra- 
diative transfer into the simulations using approximate 
schemes (Abel, Norman & Madau 1999; Gnedin & Abel 
2001; Nakamoto, Umemura & Susa 2001; Ciardi, Stoehr & 
White 2003), these simulations have emphasized the prop- 
agation of H I ionization fronts and the predicted mean op- 
tical depths. Except for the simulations of Gnedin & Abel, 
they have not solved fully self-consistently for the combined 
gas dynamics and radiative transfer. 

In this Letter, we extend the photon-conserving scheme 
of Abel et al. to include the treatment of helium and ex- 
plore the thermal effects of radiative transfer in the inhomo- 
geneous medium predicted in a ACDM cosmology. 



2 RADIATIVE TRANSFER 
2.1 Simulation data 

We use the data from a previously run ACDM simula- 
tion used to investigate effects of radiative transfer on the 
Lya forest (Meiksin & White 2003b). The simulation was 
run using a pure Particle Mesh (PM) dark matter code, and 
it was assumed the gas and dark matter have the same spa- 
tial distribution. A description of the parallel PM code is 
given in Meiksin & White (2003a). The parameters used for 
the simulation are flu = 0.30, flA ~ 0.70, fit ~ 0.045, 
h = Ho/100kms~^ Mpc"^ = 0.70, and slope of the pri- 
mordial density perturbation power spectrum n — 1.05. 
The model is consistent with existing large-scale structure, 
Lya forest flux distribution, cluster abundance and WMAP 
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where we have divided the frequencies into 100 discrete 
groups g evenly spaced between v™ and lOi'i'. (While we 
found this spacing to be adequate, we have made no attempt 
to optimize it.) These arc used in the following set of cou- 
pled equations to solve for the positions of the three I-fronts 
over time: 

driHii 



dt 



nHirm — nenHiiaHiCT) 



O-TlHII, 

a 



(5) 
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Figure 1. The baryon density fluctuations at z = 6 against co- 
moving distance from the ionizing source for the line of sight 
considered in our radiative transfer simulations. 



constraints (Meiksin & White 2003b). The simulation was 
run using 512^ particles and a 1024^ force mesh, in a cubic 
box with (comoving) side length 25 ft~^Mpc, adequate for 
obtaining converged estimates of the Lya pixel flux distri- 
bution and flux power spectrum (Meilcsin & White 2003a,b). 
A typical line-of-sight density run at « = 6, used for the ra^ 
diative transfer calculations here, is shown in Fig. 1. 



2.2 Equations of radiative transfer 

To consider the effect of radiative transfer when modelling 
the propagation of ionization fronts (I-fronts) into the IGM, 
we extend the photon-conserving algorithm of Abel ct al. to 
include helium. The advantage of this scheme is that energy 
is conserved independently of numerical resolution, ensuring 
that I-fronts propagate at the correct speeds in our simula- 
tions. In practice, this means that larger step sizes may be 
taken on the simulation space grid without the associated 
loss of accuracy which would occur when solving the radia- 
tive transfer equation through direct numerical integration. 

The probabilities for the absorption of an ionizing pho- 
ton by H I , He I and He ii , respectively, are 



driHell 

dt 



nHelTHel -|- nenHelliaHell (T) — TlHelirHell 



pHI 

pHel 

abs 

pHell 
^abs 



PHigHeigHell [l - exp (-T^"*)] /D, (1) 
gmPHeigHeii [l " exp (-T^"*)] /D, (2) 
gmgHelPHell [l - cxp (-1"^°*)] /D. (3) 



Here we have defined the auxiliary absorption and transmis- 
sion probabilities pi = 1 — exp(— r^), qi = exp(— r^), i denot- 
ing the species being referred to, ri°* = r^' -|- t^^^ + t^^^^, 
where r* is the optical depth for a given species, and 
D = pmqHeigHeii + qmPHciqucU + <?HiqHciPHcii. These prob- 
abilities can be used to calculate the ionization rate for a 
given species per unit volume, UiVi, as follows. If in a time 
5t, StNl~^ photons enter grid zone I from zone l — l, then the 
number of photons that will be absorbed in zone / by species 
i is StNl~^ P^^g, where P^bs is the absorption probability for 
species i within zone I. The ionization rate of species i in 
zone I of volume V' is then 



+ Td^ _3«(|r+!hT) 
n dt a \S n ) 



(8) 



— nenHeliaHel(T) — 3-nHeII, (6) 

a 

rfftHelll _ jjjj^jjPjj^jj _ nenHeIIlQ:HeIl(T') — 3-nHeIII, (7) 

dt a 

where iii denotes number density, F; (s^^) is the photoion- 
ization rate per atom, Ui (T) is the total radiative recombina- 
tion coefficient, and a is the cosmological expansion factor. 
The time evolution of the gas temperature T is given 

by: 

dT _ 2(G - L) . T dn, 
dt 3kn 

where n — nm + nud + nmi + niidi + niiciii + "c, "c = 
TiHii + MHeii + 2nHciii, k is the Boltzmann constant, and 
G and L (Jm~'^s~^) are the atomic heating and cooling 
rates, respectively. The last term is the adiabatic cooling 
term resulting from cosmological expansion, which will dom- 
inate the thermal effects of gas motions at the densities 
we consider. While computing the photoionization, we as- 
sume the gas over-density (not the density) stays frozen, 
which is a good approximation on the scales relevant to the 
Lya forest (Zhang et al. 1998). As in similar previous stud- 
ios, we do not include the effects of adiabatic compression 
heating as this will only afi^ect the temperature for the rela- 
tively rare virialised halos with circular velocities exceeding 
~ 20kms-^ (Meiksin 1994; Meiksin & White 2003b). 

The heating rate G' in cell I is due to the photoioniza- 
tion of H I , He I and He ii according to — Gin + Ghci + 
Ghcii, where for each species i, G\ is evaluated in a similar 
manner to the ionization rate per unit volume as given in 
Eq. (4), 



(9) 



where Xi is the ionization potential of species i. The atomic 
cooling rate L includes radiative recombination cooling to 
H I , He I and He ii , electron excitation of H i , and Comp- 
ton cooling oil the Cosmic Microwave Background photons. 
The radiative recombination and cooling rates are taken 
from Meiksin (1994). 

Eqs. (5)-(8) are solved using explicit forward Euler in- 
tegration. Although an implicit numerical scheme such as 
backward Euler integration results in a more stable solution 
at low numerical resolution (Anninos et al. , 1997), we find 
that to obtain the required numerical accuracy as well as to 
maintain stability, both methods require similar numerical 
resolution. For the sake of speed and simplicity, the forward 
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method is preferred. The timestep is restricted to being no 
more than several times greater that the hydrogen ionization 
timescale, T^j , for numerical accuracy. 

We take the frequency specific luminosity of the ionizing 
source in our simulations to bo described by a power law 
spectrum with index a = 1.5, and a luminosity of = 
10^^ WHz^^ at the Lyman edge. This spectrum is typical 
of that of QSOs, which are probable candidate sources for 
contributing to the re- ionization of the IGM. The assumed 
mass fraction of helium in the IGM is F ~ 0.235. 

We tested the photon-conserving algorithm on the pho- 
toionization of gas with uniform density around a point 
source. The resulting solutions for the ionization and tem- 
perature profiles were accurate to better than 10 per cent 
for optical depths at the Lyman edge of up to Ati, ~ 20 
per cell on the space grid. This results in a significant re- 
duction in the computational resources compared with a 
radiative transfer scheme in which the ionization rates are 
computed by direct numerical integration over the inten- 
sity (Ljje^^'^ /[47r'r]'^) and photoionization cross-sections (eg, 
Madau, Meiksin & Rees 1997). Typically, we find that to 
obtain an accuracy comparable to the photon-conserving 
scheme at /^t™ ~ 20, a direct integration algorithm re- 
quires l^T™ ~ 1/4 in each cell. For the simulation results 
presented in this Letter, the original PM grid of 1024 sepa- 
rate cells was used over a line of sight 25h~^ comoving Mpc 
in length. This provided accurate spatial resolution, as at 
most At™ ~ 6.8 in each cell at the average baryonic den- 
sity assumed for the neutral IGM at 2: = 6. 



3 RESULTS 

We set up a problem of a QSO source of luminosity at 
z = 6 placed at a comoving distance of lOft"'^ Mpc from the 
left edge of the density run shown in Fig. 1, and assume the 
gas surrounding the QSO up to that point has already been 
ionized. Placing the source at this distance avoids having to 
impose the restriction that the I-fronts propagate no faster 
than the speed of light. We indicate the displaced position of 
the source in the figures by starting the (comoving) spatial 
axes at R = 10h~^ Mpc. 

Two cases were compared: computing the ionization 
and temperature profiles by incorporating radiative trans- 
fer as discussed in section 2.2 (hereafter the RT simulation), 
and a second case neglecting radiative transfer, assuming 
instead an instantaneous uniform ionization rate across the 
whole line of sight (hereafter the nRT simulation), as is usu- 
ally done in simulations of the Lya forest. Our intention is to 
compare the final gas temperatures to determine how large 
an effect including radiative transfer may have. 

Fig. 2 shows a comparison between the IGM tempera- 
tures at z = 5 computed with and without radiative trans- 
fer. The inclusion of radiative transfer results in a signifi- 
cant boost to the temperature of the ionized IGM, a point 
illustrated by Abel & Haehnelt (1999) in the context of a 
smooth IGM. The positions of the H 11 and He 11 I-fronts at 
« = 5 in the RT simulation are 24h~^ comoving Mpc from 
the source, while the Ho iii front lags behind at about 19/i~^ 
comoving Mpc; their speeds of propagation are limited by 
the reduction in the ionization rate due to the increasing 
optical depth through the box. Consequently, the heating of 
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Figure 2. Comparison of the IGM temperature at z = -5 against 
comoving distance from tlie ionizing source. The solid line is com- 
puted using radiative transfer, wliile tlic dashed line result as- 
sumes a uniform ionization field through the entire line of sight. 




Figure 3. Comparison of the IGM temperature at 2 = 3 against 
comoving distance from the ionizing source. The line types are as 
in Fig. 2. 

the IGM, which we find to be dominated by the photoion- 
ization of neutral hydrogen at the H 11 I-front, is also con- 
strained to lie behind the H 11 I-front. All gas beyond 24/i~^ 
Mpc is still much cooler. In contrast, the temperature com- 
puted by the nRT simulation is much lower and the heating 
extends across the whole line of sight. The temperature is 
also much less sensitive to variations in the gas density for 
the nRT simulation. Fig. 3 compares the results obtained at 
« = 3. By this time, both the H 11 and He 11 I-fronts have 
reached 35/i~^ Mpc from the source. The material which has 
been ionized early on in the RT simulation is now beginning 
to cool to the same temperature as that computed by the 
nRT simulation, although the temperature differences re- 
main large. In particular, the gas that was pliotoionized at 
z = b {R < 24/1"^ Mpc), remains nearly twice as hot as 
the uniformly ionized gas by z = 3. By contrast, the regions 
which were ionized more recently are much hotter than for 
the nRT sinmlation. We find that only by z = 1, once the 
gas has been in an ionized state for a sufficiently long period 
of time, do the temperatures predicted in the two scenarios 
converge. 
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Figure 4. Plot of the energy per ionization for HI (solid line) 
Hcl (dashed line) and Hell (dot-dashed line) at z = 5. The bold 
lines correspond to the results computed using radiative transfer, 
while the lighter (flat) lines correspond to the uniform ionization 
field. 



The reason for the temperature boost in the RT sim- 
ulation is made clear by considering the energy per ioniza- 
tion {E = G/T) at the I-front. In Fig. 4, we show that 
the average energy of a photon absorbed by either H i , 
He I or He ii increases at the relevant I-front, as proportion- 
ally more of the lower energy photons have been absorbed 
already - a consequence of including radiative transfer. This 
increase to the heating energy per ionization is responsible 
for the temperature boost. The assumption of an instanta- 
neous ionization field across the line of sight results in an 
underestimate of the IGM temperature at high redshift. 

A comparison of the IGM temperature-density rela- 
tion computed by the two different simulations is made in 
Figs. 5-7. Only data points where He ii has been 90 per 
cent ionized are plotted, so that gas which is neutral or in 
the process of being ionized is not considered. In Fig. 5, 
the underdense regions in the line of sight (Fig. 1), axe ini- 
tially heated to a higher temperature in the RT simulation. 
A non-monotonic temperature-density relation results. The 
relation actually splits into two trends. To understand the 
origin of the split, we divide the data points into two dis- 
tinct sets corresponding to two regions in which the average 
density fluctuations are slightly higher or lower than the 
overall average in the line of sight (see caption to Fig. 5). 
The gas in the denser region, which is closer to the source 
and thus is ionized sooner, shows the same overall trend as 
the gas in the less dense region but is slightly cooler at a 
fixed over-density. In contrast, in the nRT simulation the 
temperature-density relation is both monotonic and essen- 
tially single- valued. By z = 3 (Fig. 6), the ionized gas has 
cooled adiabatically, with the overdense regions cooling at 
a slower rate than the underdense regions, as the denser re- 
gions are better able to maintain thermal balance because 
of the shorter cooling (and photoionization heating) times. 
The IGM temperatures computed by the two simulations are 
now beginning to converge at high density, but at low den- 
sity there is still a significant difference between the results. 
A third set of data points is also evident, corresponding to 
the latest locally underdense region to be ionized. 

Finally, by ^; = 1, the gas temperatures computed by 
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Figure 5. Comparison of the IGM temperature-density relation 
at z = 5 for regions with fully ionized helium (n£re////"-ife > 
0.9). The triangles correspond to the uniform ionization case. The 
data points for the results incorporating radiative transfer are 
split into two sets: -|- corresponds to gas in the range 
14/i~^, for which the mean over-density (p/p) = 1.452, (where p 
is the global mean gas density); and * designates gas in the range 
14/1-1 < R< 23.6/1-1, with mean over-density 0.519. 
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Figure 6. Comparison of the IGM temperature-density relation 
at z = 3. The symbols are as in Fig. 5. We also add a third point 
X designating gas in the range 23.6h~^ < R < 25.8h~^, with 
mean over-density 0.671. 

the RT and nRT simulations have nearly converged. The gra- 
dient of the line in the temperature-density plane indicates 
that adiabatic cooling due to expansion is largely responsi- 
ble for cooling the gas. We find the entropy is now nearly 
constant throughout the line of sight. The striations in the 
RT data points have also disappeared. By z = 1, the results 
from the two simulations are in good agreement. 



4 CONCLUSIONS 

We extend the probabilistic method of radiative transfer de- 
scribed for hydrogen by Abel et al. (1999) to include helium 
and tested the method on the propagation of an ionization 
front through a uniform medium. We compare the method 
against direct numerical integration of the exact time-steady 
radiative transfer equation. We confirm their conclusion that 



© 0000 RAS, MNRAS 000, 000-000 



Radiative transfer through the Intergalactic Medium 5 




iro r , , , , I , , , , I , , , , I , , , 
^— 2—1 O 1 3 

log (^/<P>) 

Figure 7. Coiriparisoii of the IGM temperature-density relation 
at z = 1. The symbols are as in Figs. 5 and 6. 



the probabilistic method accurately reproduces the position 
of the ionization fronts even for large optical depths at the 
photoelectric edges. Specifically, we recover the numerically 
integrated converged solutions for the ionization fractions 
and temperature to an accuracy of 10 per cent for an incre- 
mental optical depth per grid zone at the H i photoelectric 
edge as high as 20, while direct numerical integration of the 
radiative transfer equations requires an incremental optical 
depth less than 1/4 to obtain comparable accuracy. 

We use the probabilistic method to compute the prop- 
agation of an I-front generated by a QSO source through 
the IGM, using the density distribution drawn from a nu- 
merical simulation in a ACDM universe. Including the ef- 
fects of radiative transfer results in a substantial boost in 
the post-ionized gas temperature compared with the case of 
sudden uniform ionization with no radiative transfer. The 
boost is due chiefly to the increased energy per pliotoioniza- 
tion within the ionization fronts, a consequence of the higher 
probability for high energy photons to pass through the gas 
compared with lower energy photons. 

A consequence of the radiative transfer is to alter the de- 
pendence of temperature on density from the approximately 
polytropic relation found for moderate to low densities as- 
suming uniform photoionization (Meiksin 1994; Gnedin & 
Hui 1998) to a generally non-monotonic trend. While the dif- 
ferences between the temperatures for over-densities above 
^ 5 are small (less than 25 per cent), at lower densities the 
two temperatures converge only slowly with time, nearly re- 
covering the polytropic relation only by 2: « 1. 

At z > 1, the relation between temperature and den- 
sity is not single-valued: multiple temperatures may result 
for different gas parcels of the same over-density. Although 
this is in part a consequence of the delay in time for gas 
at different distances from the source to be photoionized as 
the ionization front propagates, an additional contributing 
factor is the environment of the gas parcel: if the ionization 
front has passed through a region of dense gas, sufficient to 
drive the optical depth at the H i photoelectric edge to near 
unity, before reaching the parcel, the radiation field will de- 
liver a higher energy per photoionization to the parcel. Thus 
the determination of the gas temperature is strongly non- 
local due to the dependence of the heating rate on the optical 



depth to photoionizing photons between the gas parcel and 
the source of ionizing radiation. 

Although the inclusion of radiative transfer results in 
a substantial and persistent boost in gas temperature, it 
is insufficient in itself to account for the broader mea- 
sured Lya absorption systems than predicted by numeri- 
cal simulations assuming sudden uniform photoionization. 
Meiksin, Bryan & Machacek (2001) find that a boost of 
AT « 1.7x 10^ K is required at z ^ 3.5, larger than the tem- 
peratures we find here. Although it may be possible to ob- 
tain a larger boost for an alternative re-ionization scenario, 
such a large boost may still require invoking late He 11 re- 
ionization. Our principal conclusion here is that for numeri- 
cal simulations to accurately predict the temperature of the 
IGM, they umst include the effects of radiative transfer. We 
are currently extending our methods to incorporate radiative 
transfer with ray tracing (eg, Abel & Wandelt 2002), into 
fully self-consistent simulations to explore this and related 
issues, including alternative sources of re-ionization like low 
luminosity AGNs and stars. 

A.M. thanks the University of Edinburgh Development 
Trust for its financial support. 
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